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Several recent experiments have established by measuring the Mandel Q parameter that the 
number of Rydberg excitations in ultracold gases exhibits sub-Poissonian statistics. This effect is 
attributed to the Rydberg blockade that occurs due to the strong interatomic interactions between 
highly-excited atoms. Because of this blockade effect, the system can end up in a state in which 
all particles are either excited or blocked: a jamming limit. We analyze appropriately constructed 
random-graph models that capture the blockade effect, and derive formulae for the mean and vari¬ 
ance of the number of Rydberg excitations in jamming limits. This yields an explicit relationship 
between the Mandel Q parameter and the blockade effect, and comparison to measurement data 
shows strong agreement between theory and experiment. 
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Ultracold gases with atoms in highly excited states 
have attracted substantial interest over recent years, for 
example for their potential application in quantum com¬ 
puting [1-3], and for study of non-equilibrium phase tran¬ 
sitions [4]. These atomic systems exhibit complicated 
spatial behavior due to strong van der Waals or dipolar 
interactions between neighboring atoms, which has been 
demonstrated through several experimental observations 
of reduced fluctuation in the number of excitations in 
ultracold gases of Rydberg atoms [5-10]. 

In these experiments, a laser facilitates excitation of 
ultracold atoms into a Rydberg state. After some time t , 
information on the mean and variance of the number of 
excited particles X(t) is obtained by repeating counting 
experiments, and the Mandel Q parameter [11] 


Q(t) 


Var [X(t)] 

~WW 


( 1 ) 


is calculated to quantify a deviation from Poisson statis¬ 
tics, since if X(t) is Poisson distributed, Q(t) = 0. The 
experiments establish that X{t) is underdispersed, i.e. 
Q(t) < 0 for t > 0, and X(t) is said to have a sub-Poisson 
distribution. 

The observed negative Mandel Q parameter is at¬ 
tributed to the Rydberg blockade effect [1, 2], There ex¬ 
ist simulation techniques [12] and models based on Dicke 
states [6] that numerically describe the Mandel Q param¬ 
eter, but to the best of our knowledge, no closed-form 
expression is available in the literature that describes the 
relation between the Mandel Q parameter and the block¬ 
ade effect. 

Explicit formulae for the Mandel Q parameter are dif¬ 
ficult to obtain, because the problem at hand is reminis¬ 
cent of continuum random sequential adsorption prob¬ 
lems [13]. The standard two-dimensional continuum ran¬ 
dom sequential adsorption problem is that of throwing 
disks of radius r > 0 one by one randomly in a two- 
dimensional box, such that the disks do not overlap. This 
process continues until a jammed state is reached. Except 


for the one-dimensional variant, such problems are noto¬ 
riously challenging to analyze due to spatial correlations. 
One further question is whether such stochastic processes 
are suited to explain effects occurring in ultracold Ryd¬ 
berg gases, and if so, under what conditions. This matter 
is discussed in [14], where a suitable stochastic process 
is provided based on rate equations that adequately de¬ 
scribe the Rydberg gas when an incoherent process (such 
as spontaneous emission) occurs [15]. 

This Letter adopts the stochastic process in [14] that 
models the Rydberg gas, and uses it to study the Man- 
del Q parameter in the jamming limit which occurs 
when atoms only transition from the ground state to 
the Rydberg state. The model includes the blockade ef¬ 
fect through so-called interference graphs, and by con¬ 
sidering specially constructed large Erdos-Renyi (ER) 
random graphs [16] that retain essential features of the 
blockade effect, we overcome the mathematical difficul¬ 
ties normally involved with having a spatial component. 
The problem remains nontrivial though, and we point 
interested readers to our rigorous derivation of the nec¬ 
essary fluid and diffusion limits [17]. This Letter explains 
how to use our theoretical insights in the context of Ryd¬ 
berg gases through less complicated heuristic arguments, 
and while doing so explicitly relates the mean and vari¬ 
ance of the number of excitations to the blockade effect. 

We consider a gas of ultracold atoms in an excitation 
volume ycR 3 , and we assume that each particle has its 
own distinct position. Each particle can go from a ground 
state to a Rydberg state, and a particle in the Rydberg 
state prevents neighboring particles from also entering 
the Rydberg state. The density of particles is assumed to 
be p , and the number of excitable particles N within any 
region A C V to be Poisson distributed with parameter 
pA. This implies in particular that in the absence of 
blockade effects, the number of excited particles within 
the excitation volume, X(t), will be Poisson distributed, 
as is the case in experiments [5-10]. It also implies that 
the particles are uniformly distributed at random over 
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the excitation volume. 

The blockade effect will be modelled using the notion 
of a blockade radius r. This is in line with simulations 
and measurements of pair correlation functions between 
atoms in the Rydberg state, which show a sharp cutoff 
when plotted as a function of the distance between the 
atoms [18, 19]. Particles within a distance r > 0 are 
considered neighbors of each other, and neighbors each 
block the other if excited. We denote the number of 
neighbors of a particle i = 0,1,..., N within its blocking 
volume Vb.i by Bi. Asa consequence of our assumptions, 
the number of neighbors of particle i is also Poisson dis¬ 
tributed. Specifically, 


P [Bi = b] 


(pM b,*) b e-^ 

6 ! 


6 = 0 , 1 ,..., 


( 2 ) 


if Vb ,i is fully contained within V. 

We will study the number of excitations by examin¬ 
ing the asymptotic behavior of large ER random graphs. 
Each vertex of such a graph will represent one particle, so 
the set of vertices is given by V = {1,..., N}. We draw 
an edge between two particles i and j if we consider par¬ 
ticles i and j to be neighbors (particles that would block 
one another). One can construct an ER random graph by 
considering every pair of vertices (i,j) once, and drawing 
the edge between i and j with probability p , independent 
from all other edges. In order to deduce information on 
X{t) through examining the ER random graph, we need 
to match the ER random graph model to the physical 
system, and we will do so by counting and matching the 
number of neighbors. Matching the models has to be 
done via the number of neighbors, because there is no 
such notion as a physical position of a particle in an ER 
random graph. This principle, in fact, makes this math¬ 
ematical model tractable. 

The number of neighbors T?er ,7 of a particle i in the 
ER random graph is binomially distributed, I?er ~ 
Bin(7V — 1 ,p), so that for 6 = 0,1,..., iV — 1 , P[T?er i = 

b] = rrVd - p and e iberA = (n- i ) P . 

When setting p = c/N where c is some constant, we see 
that as N —> oo, the distribution converges to a Poisson 
distribution, 


lim P[_B E r,i = 6 ] = —r .—, 6 = 0 , 1 ,.... (3) 

N—> oo 01 

Comparing (3) to (2), we note that the limiting distri¬ 
bution is the same if the average number of neighbors in 
the ER random graph, c, is related to the density and 
blockade volume as c = pV t,. By setting c = pV b, we 
ensure that the particles in the ER random graph have 
the same distribution of number of neighbors as in the 
spatial problem when the number of particles N —> oo. 
Figure 1 summarizes our construction. 

Let us now describe the dynamics, illustrated in Fig¬ 
ure 2 . At time To = 0 the laser is activated, and from that 




Figure 1: (left) A spatial Poisson point process in which neigh¬ 
bors within radius r block each other is used to choose appro¬ 
priate parameters for (right) an ER random graph so that the 
particles have the same distribution of number of neighbors 
as N —► oo. 


point onward excitations can occur. At a time Tj > To, 
the first particle (say 1) excites and enters the Rydberg 
state. Due to the Rydberg blockade, particle 1 will sub¬ 
sequently prevent all other particles within a radius r 
from also exciting. Later, at a time T 2 > Tj, a second 
particle excites (say 2 ), which cannot be within distance 
r of particle 1. Particle 2 from that point onward also 
blocks particles within a distance r of itself. This process 
continues until some finite time Tx(oc) < 00 when all par¬ 
ticles are either blocked or excited. The random number 
of excited particles 1 < X(oo) < N is then detected. 



Figure 2: (left) A random first particle excites, (middle) Sub¬ 
sequently, random second and third particles excite, (right) 
The process continues until all particles are either blocked or 
excited, and the resulting state is a jamming limit. 


We will now derive expressions for the Mandel Q pa¬ 
rameter. Let X m denote the excited particles at time 
T m , and let U m be the unaffected particles at time T m . 
At time T 0 = 0, no vertices are excited or blocked, so 
Xo = 0 and Uq = V. At time T m+ 1 , a uniformly randomly 
chosen particle v m +i £ U m excites, and starts block¬ 
ing random neighbors that were thus far unaffected, say 
^ 771 + 1 , 1 , - ■ •, , so that A m _|_i — X m U } , 

and — ht rn \ (} U {'Um+i.i, ■ • ■, ^m+ 1,6 })• This 

stochastic process continues until the moment r a jam¬ 
ming limit occurs, i.e. when U T = 0 and the number of 
unaffected particles equals zero. 

The number of unaffected particles U m = \U m \ can be 
described using a stochastic recursion. When the (?n+l)- 
th excitation occurs, the number of unaffected particles 
decreases by (i) the one particle that excites, and (ii) a 
random number of unaffected particles that each is neigh- 
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bor of the new excitation with probability p and thus now 
become blocked. Conditional on there being N = n par¬ 
ticles in the excitation volume, we have 

U m +i = U m -l-Bm(U m -l,p), U 0 = n. (4) 

We will now analyze the stochastic recursion in (4), and 
identify the moment r the number of unaffected particles 
is zero, i.e. U T = 0. Precisely at this moment, we have 
that the number of excitations X(oo) = r. 

Our supplemental material details the following steps 
[20]. From (4), we obtain a closed-form expression for 
E[[/ m ] by invoking the tower property and giving an in¬ 
duction argument. Through decomposition, we subse¬ 
quently obtain an expression for Var[J7 m ]. When scaling 
the probability of being neighbors as p = c/n, the mean 
and variance converge to fluid limits, which can be seen 
by letting / £ [0,1], and proving that as n —» oo, 


2.5% for the variance, and (iii) 0.015% for the Mandel 
Q parameter. Because the mean and variance are both 
overestimated, the Mandel Q parameter happens to be 
more accurately approximated. The errors our approx¬ 
imation makes can be attributed to the fact that par¬ 
ticles in the Erdos-Renyi random graph model have no 
physical position, whereas particles in two-dimensional 
Poisson disk throwing processes do exhibit spatial cor¬ 
relations. Intriguingly the random graph, which has no 
spatial interpretation, yields a good approximation. 



E [C/[/„]] 


-t u{f), 


Var[f7 [/n] ] 

n 


-t v{f). 


(5) 


Here, [■] denotes rounding to the nearest integer, and the 
fluid limits are u(f) = e~ c f — (1 — e -c ^)/c, and v(f) = 
(e -c ^(l — e -c ^)((l + 2 c)e~ c f — l))/(2c). Note that these 
fluid limits are rigorously proven in [17]. 

Consider now Figure 3 (left) and the following steps. 
The process U m hits zero when m ~ f*n, with f* = 
ln(l + c)/c being the solution to u(f*) = 0. Therefore, 


E[X(oo)\N}* f*n= nln(1 + c) , (6) 

c 

To approximate the variance, calculate u'(f) and note 
that u'{f* +e) ~ — 1 for sufficiently small e. Since U m is 
probably near 0 for m ss f*n, the fluctuations in A'(oo) 
will thus be of the order of ^/Var[{/[/*„]]. Hence, 


Var[X(oo)|W] « Var^,^] « v{f*)n = + (7) 


Figure 3: (left) The fluid limit u(f) together with a sample 
path of Um/n for n = 50. The dashed line indicates the 
tangent at /*, and the arrows indicate the typical fluctua¬ 
tions. (right) Histogram of the number of excitations in a 
two-dimensional random sequential adsorption problem, and 
precisely n = 800 particles, together with the probability den¬ 
sity function of a normal distribution with mean n In (1 + c)/c 
and variance nc/ (2(1 + c) 2 ). 

It is important to understand that the results thus far 
are conditional on there being N = n particles within 
the excitation volume. However, the number of parti¬ 
cles N ~ Poi(pU) is itself random. To obtain an un¬ 
conditional expression for the mean and variance, we can 
utilize the tower property, E[X(oo)] = E[E[X(oo)|lV]] ~ 
E[iV] In (1 + c)/c = pV\n (1 + c)/c, and decomposition, 
Var[X(oo)] = E[Var[X(oo)|IV]] + Var[E[A(oo)|fV]] « 
E[fVc/(2(1 +c) 2 )] + VarfiV In (1 + c)/c] = (c/(2(l + c) 2 ) + 
(ln(l + c)/c) 2 )pV. Recalling definition (1), the Mandel 
Q parameter in the jamming limit is therefore 


Invoking the central limit theorem, we have for large fixed 
n that the number of excitations is approximately normal 
distributed with mean nln(l + c)/c and standard devi¬ 
ation Y / nc/( 2(1 + c) 2 ). This, (6), and (7) are formally 
established by deriving diffusion limits in [17]. 

Let us compare (6) and (7) to simulations of the mean 
and variance observed in the two-dimensional random se¬ 
quential adsorption problem described earlier, and with 
periodic boundary conditions. We consider h = 1pm, 
l = w = 400pm, r = 6.5pm, and p = 5 x 10 9 cm -3 , 
which are typical values in magneto-optical traps, and 
correspond to n « 800 and c ~ 0.664. Figure 3 (right) 
shows a histogram of the number of excitations, as well 
as the probability density function of a normal distribu¬ 
tion with mean nln (1 + c)/c and variance nc/(2(l + c) 2 ). 
Comparing to the simulation’s outcome, our expressions 
differ for this set of parameters (i) 2.6% for the mean, (ii) 


Q(oo) 


c 2 ln(l + c) 

2(1 + c) 2 In (1 + c) c 


( 8 ) 


which is exact in the ER case when pV —► oo [17]. Note 
that (8) only depends on the average number of neigh¬ 
bors c, which in fact explains observations on simulated 
Mandel Q parameters [12] as we discuss in [20]. 

Let us also discuss the time-dependency of the mean 
number of excitations. We incorporate time-dependency 
by assuming that every unaffected particle excites at rate 
A, and specifically that T m —T m _i ~ Exp(A17 m _i), which 
corresponds to modelling the Rydberg gas using rate 
equations [15] as discussed in [14]. Under these assump¬ 
tions, we obtain the time-dependent fluid limit [17] 


E[X(f)|7V] 

n 


v{t) = A f 
Jo 


i(x(s))ds 


(9) 
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After substituting u(f) = e~ c f — (1 — e~ c ^)/c into (9), 
recalling that initially no particles are excited, and taking 
the derivative, we obtain the differential system dx/d t = 
A(exp (— cx(t )) — (1 — exp (— cx(t)))/c), with x(0) = 0, for 
x(t). This differential system has as its unique solution 
x(t) = In (1 + c — ce _At )/c, and in particular, we recover 
the mean fraction of excitations in the jamming limit by 
calculating lim^oo x(t) = In (1 + c)/c. 

We now validate our model by comparisons with ex¬ 
perimental data in [ 6 , 7], which requires us to incorpo¬ 
rate the notion of a detector efficiency 77 £ [ 0 , 1 ] into 
the model. The detector efficiency 77 can be interpreted 
as being the probability that a Rydberg atom is de¬ 
tected. Let /, ~ Ber(? 7 ) denote random variables that 
indicate whether each i-th Rydberg atom is detected. 
The number of detected Rydberg atoms is then given 
by X D (t) = Y^i=i !i- Assuming the h,... , I X {t) are in¬ 
dependent, calculation shows that E[Wo(t)] = ? 7 E[X(f)], 
and Var[A'£>(f)] = 77 2 Var[X(i)] + 77(1 — 77 )E[X(f)], see our 
supplementary material [20]. The detected Mandel Q 
parameter thus reduces to Qo{t) = 77Q(i), see also [5]. 

The experiments in [6] were on excitation volumes said 
to contain pV = 8 x 10 3 ground-state atoms, and with a 
reported detector efficiency of ?/ = 0.40. Fitting 


E[X D (t)] * wrMl+c-oe-") 


( 10 ) 


to measurements of the number of excitations as a func¬ 
tion of time [ 6 , Fig. 1(a)], we obtain an excitation rate 
of A = 14kHz, and average number of neighbors of 
c = 2.7 x 10 2 . Figure 4 shows strong agreement between 
theory and experiment. 



Figure 4: The average number of detected excitations as a 
function of time, E[A£>(t)], fitted to the measurement data 
in [6, Fig. 1(a)]. The fit results in an excitation rate of A = 
14kHz, and average number of neighbors of c = 2.7 x 10 2 . 

Lastly, we will compare our model to a histogram of the 
number of detected dark-state polaritons in [7]. The his¬ 
togram displays sub-Poissonian statistics due to a block¬ 
ade effect that is a result of the dominant Rydberg char¬ 
acter of the polaritons. Because of a partial overlap be¬ 
tween the excitation laser and the cigar-shaped atomic 


cloud, we will infer the size of the excitation volume us¬ 
ing the density p = 5 x 10 17 m -3 [7] as follows. The 
detector efficiency is reported to be 77 = 0.4, and the 
histogram has a sample mean of E[Xd(oo)] « 11. If 
we assume that the blockade regions are spherical, and 
since the blockade radius r « 5pm [7], we find that 
c = |p 7 rr 3 « 2.6 x 10 2 . Using our formula for the 
mean number of detected Rydberg atoms, it follows that 
V = cE[A D (oo)]/(p 77 ln(l + c)) « 2.6 x 10- 15 m 3 . The 
factor with which the density function of the Poisson 
distribution is scaled in [7, Fig. 4(a)] is n s « 315. Fig¬ 
ure 5 now compares the appropriately scaled probability 
density function of a normal distribution with mean and 
variance as predicted by the model to the histogram in 
[7, Fig. 4(a)]. Our result Qd ~ —0.36 is consistent with 
their observation that Qd = —0.32 ± 0.04 in the density 
range 2 x 10 17 m -3 < p < 2 x 10 18 m -3 . 



Figure 5: Histogram [7, Fig. 4(a)] of the number of detected 
Rydberg atoms, together with the appropriately scaled prob¬ 
ability density function of a normal distribution with mean 
E[A'_d(oo)] and variance Var[A_o ( 00 )]. Here, Qd( 00 ) ~ —0.36, 
and the dashed line indicates the Poisson distribution. 

This Letter derived closed-form expressions for the 
Mandel Q parameter in limiting large random graphs 
constructed to model the spatial problem. This approach 
allowed us to derive explicit formulae for the mean and 
variance of the number of Rydberg excitations in the jam¬ 
ming limit, that turn out to be functions only of the av¬ 
erage number of neighbors within the blockade volume. 
Our comparison to measurement data of [ 6 , 7] shows 
quantitative agreement between theory and experiment, 
and we conclude that the model captures blockade effects 
observed in ultracold Rydberg gases. 

Interesting future research would be to explore the ap¬ 
proximating relation between random graphs and spatial 
problems, particularly because higher-dimensional con¬ 
tinuum random sequential adsorption processes are diffi¬ 
cult to analyse. The underlying stochastic recursions can 
also be generalized to incorporate additional effects [17], 
such as a slower reduction in the number of neighbors as 
more particles become blocked, and this can potentially 
extend the use of random graphs as an approximation to 
particle systems that exhibit complicated interactions. 
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SUPPLEMENTARY MATERIAL 


Solving the stochastic recursions 


Our exposition led to the following stochastic recur¬ 
sions for the number of such particles, 


Xm -11 — 7^1 4 “ 

U m+1 = U m -l-Brn(U m -l,p), (11) 


with Xq = 0 and Uq = n. These recursions can be lever¬ 
aged to determine the mean and the variance of the num¬ 
ber of unaffected particles at each moment m excitations 
have occurred, i.e. E[[/ m ] and Var[f/ m ]. To see this, let 
in > 1 and start by noting that 


U m = Bin(f/ m _! — 1,1 — p). (12) 

Utilizing the tower property, we find that 

E[U m ] = E[E[Bin([/ m _i - 1,1 - p)\U m ]] 

= (l-p)E[U m - 1 ]-(l-p). (13) 

Iterating and recalling that E[C/q] = n , we obtain 


E[d m ] = (l- P rn^(Wp) ! 


i—1 


= (1 -p) m n- 


(i - P ) ^ (i - P r +i 


(14) 


A recursion for the variance can be found in a similar 
fashion by first decomposing 


Var[U m ] = Var[E[Bin(f7 m _i - 1,1 - p)\U m -i]] 

+ E[Var[Bin([/ m _i - 1,1 - p)|E/ m _i]] 

= Var[(C/ m _i - 1)(1 - p)} + E[(t/ m _i - l)p{\ - p)] 

= (1 -p) 2 Var[C/ m _i] + p(l - p)(E[U m -i] - 1), (15) 
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and then recalling that Var[Z7o] = 0 , E[f7o] = n, so that Substituting (14) into (16) and simplifying, we find that 
after iterating, 

Var[U m ] = p(l — p) 2m ~ l {n — 1) 

m— 1 

+ P j2^~p) 2l -\num-i\-i). (i6) 

2 = 1 


Var [U m \ = 


(p - 2)(1 - p) m ((n - 1 )p + 1 ) + (1 - p) 2m (l - (n - l)(p - 2 )p) -p+1 

(p ~ 2 )p 


(17) 


r 


Determining fluid limits 

When n —> oo and p = c/n , there exist fluid limits for 
E[J7 m ] and Var[f7 m ]. To see this, define 


u(f) = lira E[U [fn] \/n 


(18) 


for / € [ 0 , 1 ], with [•] denoting rounding to the nearest 
integer. Utilizing (14), we find that 


, , / c\ I/”! 1/ c\ / / 1\ i/ ra J\ 

i(f) = lim (l-)-(l-) (l — (l-) ) 

n—too \ n/ C\ 71/ \ \ C/ / 


= e" c/ - -(l-e- c/ ). 
c 

Similarly, after defining 

v(f) = lim Var[[/[/„]]/n, 

n—>oo 1 

and substituting (17), we find that 

e" c/ (l - e" c/ )((l + 2 c)e~ cf - 1 ) 


1 \ [/™] > 
Cj 

(19) 


( 20 ) 


v(f) = 


2 c 


( 21 ) 


Incorporating detector efficiency 

Assuming the ,Ix(t) are independent, the aver¬ 

age number of detected Rydberg atoms is 


x(t) 

E[X D (t)]=E[Y,Ii 


i=1 


= E 


x (t) 


E[^A|A'(t)]] =77E[X(t)]. (22) 

i=1 


The variance of the number of detected Rydberg atoms 
is 

x(t) 

Var[X D (t)] = Var R 


i =1 


= Var 


x(t) 


E 




-E 


x(t) 

Var[]T Ii\X(t) 


X{t) 


ry 2 Var [X{t)] + E Varft] 


2=1 


= r/ 2 Var[X(t)} + 77(1 - ry)E[X(f)]. 

The detected Mandel Q parameter is therefore 

which completes the derivation. 


(23) 


(24) 


Comparison to Petrosyan, 2013 

Reference [12] describes usage of semiclassical Monte 
Carlo simulations to study stationary states of the Ryd¬ 
berg gas in a two-dimensional system. There, parti¬ 
cles are positioned on points of a lattice with spacing 
a = 532nm, that fall within a circular excitation area 
with radius R , which is varied relative to the blockade 
radius of r ss 1.905pm. Ref. [12] finds numerically that 
Q ss —0.84 for R > r. This independence on the sys¬ 
tem size is explained by our model, because ( 6 ) and (7) 
indicate that the Mandel Q parameter only depends on 
the average number of neighbors c, which in this simula¬ 
tion setup approaches a constant for sufficiently large R. 
Modelling the blockade area as a hard circle of radius r, 
we have by Gauss’s circle problem that c+ 1 ss 37 lattice 
points fall within the blockade area for sufficiently large 
R. Because the number of particles within the excitation 
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area N did not fluctuate between simulation instances, suits in Q « c 2 /(2(l + c) 2 ln(l + c)) — l | c =36 = —0.87, 

we can use the conditional expressions for the mean and which is close to the simulation result, 

variance to estimate the Mandel Q parameter. This re- 



